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We derive an improved rigorous lower bound on the space and time averaged 
temperature (T) of an infinite Prandtl number Boussinesq fluid contained be- 
tween isothermal no-slip boundaries thermally driven by uniform internal heat- 
ing. A novel approach is used wherein a singular stable stratification is intro- 
duced as a perturbation to a non-singular background profile, yielding the estimate 
(T) > 0.419 [R logR] _1//4 where R is the heat Rayleigh number. The analysis relies 
on a generalized Hardy- Rellich inequality that is proved in the appendix. 
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I. INTRODUCTION 

Rayleigh-Benard convection, the buoyancy-driven flow of a fluid heated from below and 
cooled from above, is a fundamental paradigm of complex nonlinear dynamics, pattern 
formation, and turbulence. An ongoing challenge for analysis, theory, computation, and 
experiment is to ascertain how the heat transport depends on the thermal forcing as gauged 
by a nondimensional Rayleigh number and the fluid's material properties (typically the 
Prandtl number, the ratio of kinematic viscosity to thermal diffusivity) pQ . Rigorous bounds 
on heat transport in Rayleigh-Benard convection within the Boussinesq approximation were 
pioneered by Howard [2] and subsequently elaborated by Busse [3j- Later, following the 
motivational work of Hopf [I], an alternative variational framework for bounds on turbulent 
transport of momentum, mass, and in the case of Rayleigh-Benard convection, heat, known 
as the 'background method' was formulated [3]. This is the approach we adopt here. 
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In this paper we consider an infinite Prandtl number Boussinesq fluid contained between 
two rigid isothermal boundaries thermally driven by constant internal heating. This model 
is inspired by convection in the Earth's mantle where the Prandtl number is O(10 24 ) and 
the motion is predominantly driven by a semi-uniform heating from radioactive decay. For 
defmiteness we consider an isoviscous fluid subject to no-slip isothermal vertical boundary 
conditions (without loss of generality T = on the boundaries) with uniform heating in 
the bulk, an idealization of the actual geophysical conditions. At low heating rates as 
measured by the dimensionless 'heat Rayleigh number' R, the fluid remains at rest and heat 
is transported to the boundaries by conduction with a parabolic temperature profile across 
the layer. At high R convection sets in, increasing the effective thermal conductivity of the 
layer and consequently lowering the bulk averaged temperature. The enhancement of heat 
transport is indicated by (T), the ratio of the bulk averaged temperature to that of the 
parabolic profile resulting from conduction alone. The challenge is to determine how (T) 
varies with R. The no-convection conduction solution exists (be it stable or unstable) for 
all R and realizes the upper limit on the bulk averaged temperature, (T) < 1. At high R 
the question is, how low can (T) go? 

This problem was previously considerd by Lu et al [6] who used estimates originally 
derived for boundary driven convection [7] and a simple piecewise linear background profile 
to produce a lower bound on the space and time averaged temperature relative to that of 
the conduction state. That result was 

(T) > 0.81 R~ 2/7 . (1) 

Because 2/7 = 0.2857. . . , this is not inconsistent with a scaling law measured from direct 
numerical simulations [H] suggesting that 

(T) ~ 1.65 R-°- 234 . (2) 

We note, however, that those computations employed 'free-slip' boundary conditions on 
the velocity field rather than the no-slip conditions employed in the analysis. Boundary 
conditions for traditional Rayleigh Benard convection can drastically affect the dynamics so 
the comparison must be taken with a degree of caution. 

More recent developments [9j [10] indicate that background profiles including some stable 
stratification in the bulk may be optimal for infinite Prandtl number convection suggesting 
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room for improvement for rigorous results. In particular, a singular integral analysis was 
performed to obtain a key estimate that was then utilized in the background method to 
produce an upper bound on the Nusselt number Nu, the dimensionless measure of the 
enhancement of heat transport in boundary-driven Rayleigh-Benard convection, in terms of 
the traditional Rayleigh number Ra of the form Nu < [Ralog(Ra)] 1//3 [10J. Here we show 
that that key estimate is in fact a modified Hardy-Rellich inequality and we derive the sharp 
prefactor. The newly derived inequality is then applied to the internal heating problem via 
the background method and, along with some additional considerations, we prove 

(T) > 0.419 R" 1/4 (log R)" 1/4 . (3) 

The rest of this paper is organized as follows. The next section describes the Boussinesq 
equations of motion with internal heating and provides an outline of the background method 
applied to the problem. Section III introduces the particular background temperature field 
as a logarithmic perturbation of a quadratic profile and applies the modified Hardy-Rellich 
estimate to obtain the bound (|3]). Section IV discusses these results and briefly remarks on 
the parallels between the internal heating and boundary driven convection problems. The 
new derivation of the Hardy-Rellich inequality is described in the appendix. 



II. INTERNAL HEATING AND THE BACKGROUND METHOD 

In non-dimensional variables the Boussinesq approximation of the Navier-Stokes equa- 
tions with a constant internal heat source is 

+ u-Vu \+Vp = V 2 tt + RTfc (4) 



Pr \dt 

<9T 

— + u ■ VT = V 2 T + 1 (5) 

at 

V-u = (6) 

where u = iu + jv + kw is the velocity field, T is the temperature field [BJ. The heat 
Rayleigh number is R and Pr is the Prandtl number. We consider periodic horizontal 
boundary conditions in all variables with no-slip conditions (u = 0) at the bottom {z — —1) 
and top (z = 0) of the box. Combined with the incompressibility condition (J6|, this implies 
that dw/dz = at the top and bottom boundaries as well. For definiteness we choose the 
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special case where the boundary temperature is isothermal: 

TL = _ 1 = 0=T| z=0 . (7) 

Defining the space-time average of a function f(x, y, z, t) as 

If 1 1 f Lx 1 f Ly f° 

(/> = lim - dt' — dx — dy dzf(x,y,z,t') (8) 
t^oo t J L x Jo L y J J_ x 

(assuming that the limits exist) we are interested in obtaining a lower bound on the average 

temperature (T) in terms of the Rayleigh heat number R. We focus on the infinite Pr 

limit of Q, the rigorous validity of which has recently been established [TT] . so that the 

Navier-Stokes equations become the Stokes equations 

Vp = V 2 u + RTfc. (9) 

To apply the background method, break the temperature field up into a background pro- 
file and fluctuations according to T(x,y, z, t) = t(z) + 9(x,y,z,t) where r(— 1) = r(0) = 
leaving the fluctuation 9(x,y,z,t) to satisfy homogeneous Dirichlet conditions at the top 
and bottom boundaries as well. Applying this decomposition to the equations of motion 
and considering the space time average of both the momentum and properly weighted tem- 
perature equations we arrive at the following bound on the average temperature (see [6] for 
details of the derivation): 

(T) > 2(r) - ((rf> (10) 
as long as the quadratic (in 6) functional 

H = (|W| 2 ) + (2t'w9) (11) 

is positive semidefinite among and temperature fluctuations and velocity fields satisfying 
the boundary conditions. H is quadratic in 6 because, applying the curl operator twice to 
rt9j), we see that there is a linear albeit nonlocal instantaneous slaving of w to 9 

A 2 w = -RA H 9 (12) 

where A is the full Laplacian while A# is the horizontal Laplacian in x and y. 

For calculational convenience we apply the Fourier transform in the horizontal directions 
to obtain the relation for each wave number k = |k|, 
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where now for all k the single-wavenumber quadratic forms 



dz 



+ k 2 \9 k \ 2 + 2T , Re[6 k w ]s } 



dz 



must all remain positive semidefinite. In the following we consider wavenumber by 
wavenumber so we may drop the : and subscript k. In other words we seek to maximize 
2(r) — ((V) 2 ) while maintaining positivite-semidefiniteness of 



H :-- 



dB 



dz 



k 2 \6\ 2 + 2T'Re[6w] 



dz 



(13) 



uniformly in k, where 9{z) satisfies Dirichlet boundary conditions and w(z) solves 



w"" - 2k 2 w" + k*w = Rk 2 8, 



(14) 



while satisfying both Dirichlet and Neumann boundary conditions on on [—1,0]. 

Previous analysis of this problem and boundary driven convection considered background 
profiles r(z) constant in the bulk of the layer so the only possible negative contribution to 
H relied on the product of w(z) and 6(z) in boundary layers where both are constrained to 
be relatively small in magnitude. As was shown for boundary driven convection, however, 
a stably stratified (i.e. t'(z) > 0) profile in the bulk can be exploited to utilize the positive 
weighted correlation between w and 9 resulting from the slaving and improve the positivity 
of H, allowing for sharper estimates |10j . 



III. SINGULAR PERTURBATION OF A STABLY STRATIFIED PROFILE 

Consider the family of background profiles illustrated in Figure [TJ 

, a log +6(1 - z 2 ) -\<z<-5 
r(z)={ 1 ' (15) 

=f [a log(l/5) + 6(1 - 5 2 )) -5<z<0, 

where the positive parameters 5 < 1, a, and 6 will be chosen to optimize the bound. The 
logarithmic term enhances the positivity of H, and hence leads to an improved scaling of the 
boundary layer with R, while the quadratic term is meant to increase the integral of r{z) 
sufficiently to offset the slow logarithmic growth near z — — 1 and lessening the negative 



impact of the Dirichlet integral in (10). 



6 



0.05 




FIG. 1: Background profile (15) 



It is easily verified that 



r(z)dz 



{t'{z)) 2 dz 




5 
2 



5 1 
6 



( 2 log(l/^) + 4 _ 25\og{l/5) - 46 



+ h 



(16) 

(17) 

(18) 
(19) 



producing the lower bound on the average temperature (T) given by (10) when the positivity 



of H is maintained. An appropriate choice of the scaling of the parameters a and b with 
respect to 5 will allow us to determine both the 'correct' boundary layer scaling, and to 
maximize the lower bound on (T). 

The two key inequalities required for the analysis are 



6{z)w*{z)zdz < 



-i 



and 

< 4 

i 



z 



R 



dz < 0. 



(20) 



(21) 



The first inequality (21) is an exercise in integration by parts the details of which are left 



to the reader. The second inequality (21) is a restatement — and slight improvement in the 



prefactor — of the key result previously derived for Rayleigh-Benard convection [10J. While 



7 



a prefactor improvement may be considered minor, our approach to prove it is simplified 
significantly and embeds the problem in the familiar context of a generalized Hardy-Rellich 
inequality. The proof, provided in the Appendix, also indicates that the estimate with this 
prefactor is sharp. 

To determine conditions guaranteeing the positivity of H we reformulate it neglecting 



much of the L 2 norm of j- as well as the fc 2 |#| 2 term and use (20) to observe 



H > 



-6 



dO 



dz 



dz — 2a 



Re[9w* 



-dz 



-i 



/2alog(l/5) 2a 



-s 



Abz + 



26(1 - 8 2 



Re[9w*]dz. 



(22) 
(23) 



In the above we added the bulk terms to the boundary layer in order to apply ( 20 ) and (21) 

ren implies 

/2alog(l/5) 2a 



to the entire interval. Applying (21) then implies 

, 8a \i 
dz — — I — 3 



H > 



d9 



dz 



R 



-dz 



z^ 



5 



Abz + 



z 



26(1 - 5 2 
6 



Bound the magnitude of the last integral in (|24|) as follows: 

-5 



2alog(l/<J) 2a 



5 



46^ 

z 



26(1 - 5 2 



Re[9w*]dz 



< 



/2alog(lM) 2a A1 26(1 - 5 2 ) 

1 — 7 Abz H 

z 



z 



\o\ 1 


w 




\z 


1/2 


z 


3/2 



< 2 sup 



\0(z)\ 



X 



6<z<0 \Z 




1/2 



alog(lA) a Ul-5 2 ) 

— 2bz + 

6 z 



dz 
2M 2 



dz 



\w 



2 \ 1/2 
:rdz 



\z\ 



Re[6w*]dz. 
(24) 



(25) 
(26) 

(27) 

(28) 



1/2 



The homogeneous boundary conditions on 9{z) mean that for z G (—8,0), 



\0(z)\ 



°d9 ( 

dz 



-dz 



< z 



1/2 



d9 



dz 



1/2 



dz 



< z 



1/2 





-s 



dd 



dz 



1/2 



dz 



(29) 



Hence we can bound the supremum in (27) and apply Young's inequality to see that 
r0 /2alog(l/5) 2a 26(1 - 8 2 



< 



-5 


-s 



d9 

dz 



8 



dz + 



Abz + 



8 



Re[ew*]dz 



(30) 



-<5 



alog(l/8) a 
8 z 



-2bz + 



6(1-5 



2M2 



dz x 





-1 



w 



-dz. (31) 
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Inserting (31) into (24) we see that 



H > 



R 



alog(lM) a , b(l-5 2 ) 
2fe + — '- 

Z 



dz 



\w? 



dz. 



(32) 



The integral about the boundary layer in (32) can be computed exactly At this point 



we choose a = a'6/ log(l/<5) and b = b'6 where a' and b' are 0(1) absolute constants. Then 



o 

z 4 

-S 

,5 5 



ojgg(l/g a %z | b(l-5 2 ) 
5 z ' 5 



dz 



' 5 1 
5 



,5_ 5 
5 



a - + 2a'b'- + b' 2 - + (5 5 log(l/5)) 



(33) 
(34) 



as 5 —7- 0. Comparing this with (32) we see that the minimal requirement for H to remain 



positive in the 5 — > or R — > oo limit is 



8 , £ 4 [log(l/£)] o| / 4 log(l/<S) y 2 ^ 4 log(l/(5) 

— ~ a h 2o 

R 5 5 a' 5 

1 
R 



5 4 log(l/5) 



where 



(a' + Z/) 2 
10a' 



This yields the scaling of the boundary layer thickness as 



(35) 
(36) 

(37) 



5~[£(a',&')Rlog(R)]- 



-1/4 



(38) 



The average temperature is bounded by two times (16) minus (17) implying that, asymp- 
totically, 



(T) > \b'8 - a ,2 5 - 2a'b'5 - b ,2 5 



b ' _ a « _ 2a 'b' - b' 2 ^j £{a', b')- 1/A (Rlog(R))" 



1/4 



To obtain the 'best' prefactor, we maximize over a' and b' to achieve 

o3/4cl/4 

(T) > — — (Rlog(R)) 1/4 ~ 0.419 (Rlog(R)) 1/4 
6 



(39) 
(40) 



(41) 



where the optimal prefactor is obtained for a' = jx and b' = j^. 



16" 
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FIG. 2: The background profile (15) where a and b scale in the optimal sense compared to the 



logarithmic profile given by ( 42 ) . 



IV. DISCUSSION AND CONCLUSIONS 



The background profile (15) can be considered the sum of a singular logarithmic profile 



and a smooth conduction-like quadratic profile. If the logarithmic term only is considered, 
i.e., 6 = 0, then the profile would be 



alog(^) 

alog(l/5)z 



1 < z < -S 
-5 < z < 0. 



(42) 



This would be analogous to the approach taken in [10] for boundary driven convection. 
However if the same steps are followed one sees that the optimal bound occurs for a ~ 



[i°g(i/<5)r 



in which case the bound becomes 



<T>> 



(43) 



This bound on the averaged temperature is weaker than that derived previously. However, 



the same analysis performed previously to ensure the positivity of H, but with (42) as the 
background profile yields the scaling 5 ~ R 1 / 4 . 



Fig. 2 yields some insight. The purely logarithmic profile (42), depicted as the dashed 



plot, yields a thicker boundary layer at high R because the steep gradient near z = enhances 
the positivity of H sufficiently to maintain the increased size of the boundary layer. But this 
costs dearly in the computation of the Dirichlet integral that negatively affects the estimate 
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of (T) while adding very little to the computation of (tq(z)). That is, while the quadratic 
term (b > 0) thins the boundary layer, it also contributes significantly to (r) and raises the 
lower bound. 

In previous applications of the background method [7J HOI C2] the scaling of the boundary 
layer dictates the bound: typically the heat transport is bounded by where 5 is the size of 
the boundary layer. Bounding the average temperature from below for the internal heating 
problem creates a different situation where the 'optimal' boundary layer scaling in terms of 
S yields a sub-optimal bound in terms of R. 



The bound (41) is not inconsistent with numerical simulations [S] although we reiterate 
that the simulations employed stress-free (a.k.a. free-slip) boundary conditions on the veloc- 
ity, as opposed to the no-slip conditions employed here, that may affect the scaling behavior. 
The stress-free internal heating problem is addressed in [13]. It is also of interest to consider 
numerical solutions of the Euler-Lagrange equations for this problem, as has been done for 
the boundary driven convection for finite [H] and infinite [U] Prandtl number Rayleigh- 
Benard convection. Their solution would indicate what the true optimal background profile 
is, and may provide additional insight into the pursuit of further rigorous bounds. 
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Appendix A: A generalized Hardy-Rellich inequality 



We will establish (21 ) for all functions w(z) and 6(z) that satisfy (14) with the prescribed 



boundary conditions. Note that with the change of variables z — > —z this is equivalent to 
casting the problem on the positive unit interval as 



Re 



1 6w* 4 

dz > — 

z ~ R 



1 M 2 



dz 



(Al) 



where (14) is satisfied for z e [0, 1] and w(0) = w(l) = w'(0) = w'(l) = 6(0) = 6(1) = 0. In 



this context, (Al) is recognized as a factor of two improvement on the original proof [10J. 
As in the original proof we will prove the following proposition: 
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Proposition A.l If < c < d < oo, the smooth function w(z) satisfies 

w(c) = = w(d), w'(c) = = w'(d), 



and 9(z) is defined by w"" — 2k 2 w" + k A w = Rk 2 6, then 



Re 



d 6w* 4 

dz > - 

z ~ R 



\w\ 



-dz. 



(A2) 



(A3) 



In order to see the connection between (A3) and Hardy- Rellich inequalities, make the 
change of variables w(z) = z 1 ^ 2 ((z). It follows that ((z) also satisfies (A2). Inserting this 
change of variables into the fourth order term that results from the definition of 0(z), we 
see that 



dz 



\Cfdz 



b \/-i\2 



IC 



-dz + 



45 
16 



ICP 



dz. 



A similar calculation leads to 



6 w"w* 



dz 



\C\ 2 dz + 



4 



dz. 



(A4) 



(A5) 



Putting (A4) and (A5) together, we see that (A3) can be restated as 



Lemma A.l For smooth functions ((z) satisfying the boundary conditions (A2), 



IC 



//|2 



3|CT , 45 |CP 



2 z 2 



16 z 4 



dz + k 2 



> 4k 2 



2|C1 2 -^) dz + k*J*\(\ 2 dz> 

ICI 2 

^dz. 



(A6) 



Traditionally a Hardy-Rellich inequality is formulated in terms of the LP norms of the 



operator D q = 4-^ where q = 1, 2 and possibly higher orders (see [15] for example). (A6) is, 
with the appropriate integrations by parts, nothing else than the L 2 norm of the differential 
operator D 2 — k 2 acting on ((z). The inclusion of the wave number k here causes us to refer 
to this inequality as a generalized Hardy-Rellich inequality. 

To prove the Lemma, consider the following one-parameter family of integrals, 



< 



(D 2 - k 2 ) 



c 



dz, 



(A7) 



where ((z) satisfies the homogeneous boundary conditions. Expanding (A7) and integrating 
by parts multiple times leads to the following identity: 



\C"\ 2 dz + 2v(i/-2) J ^J r dz + 2k 2 J |Cf cfe + z/(z/ + 6 + z/ 3 - 4z/ 2 ) J ^dz 



+ k 4 I \(\ 2 dz > 2v 2 k 2 



ICP 



dz. 
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Letting v = | produces 



+ k 



\C\ z dz 
a 



ICT. 45 
z A 16 



ICP 



cfe + 2k 2 



\C?dz 



9. 



ICI 2 ^ > ~k 2 



ICP 



which is easily rearranged to establish the Lemma. 

The methodology employed above is not restrictive to this particular situation, and lends 
itself immediately to extension to higher order operators, and possibly higher dimensions 
as well. The free parameter v can be adjusted as desired, indicating a significant utility 
to this method of producing Hardy-Rellich type inequalities. Hardy-Rellich inequalities 
with remainder terms can also be computed by optimizing over the wave-number k (for 
an example of other Hardy-Rellich type inequalities with remainder terms see the work of 
Evans and Lewis [TB]). 

The strictness of the inequality derived here can be verified by considering functions ((z) 



that saturate (A7), that is those functions satisfying the boundary conditions together with 



(D 2 - k 2 ) 



0. 



(A8) 



Solutions of (A8) are linear combinations of modified Bessel functions: 



CO*) 



,1/2 + !/ 



[C x K q (hz) + C 2 I q (kz)} 



(A9) 



where q = \flv 2 +~2v~+l/4. Just as the original Hardy inequality p2] is not saturated 



for any nontrivial analytic functions, the functions (A9) cannot satisfy all the boundary 



conditions simultaneously so there is no analytic solution to (A8) that saturates (A7). How 



everm regularizing (A9) appropriately at the boundaries will produce a sequence of functions 



that satisfy the boundary conditions and, in the unregularized limit, satisfy (A8). Hence 



while (A7) is never saturated, there can be no improvement on the prefactor derived by this 



method, i.e., the approach outlined here is not only robust and amenable to adaptation, but 
also produces sharp estimates. 
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